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Abstract 

A novel method to solve inverse problems for the wave equation 
is introduced. The method is a combination of the boundary control 
method and an iterative time reversal scheme, leading to adaptive 
imaging of coefficient functions of the wave equation using focusing 
waves in unknown medium. The approach is computationally effective 
since the iteration lets the medium do most of the processing of the 
data. 

The iterative time reversal scheme also gives an algorithm for con- 
structing boundary controls for which the corresponding final values 
are as close as possible to the final values of a given wave in a part of 
the domain, and as close as possible to zero elsewhere. The algorithm 
does not assume that the coefficients of the wave equation are known. 

Keywords: Inverse problems, wave equation, control, time reversal. 
AMS classification: 35R30, 93B05. 



1. Introduction 

We present a novel inversion method for the wave equation. Suppose that 
we can send waves from the boundary into an unknown body with spatially 
varying wave speed c(x) . Using a combination of the boundary control (BC) 
method and an iterative time reversal scheme, we show how to focus waves 
near a point Xq inside the medium and simultaneously recover c(xq) if c is 
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isotropic. In the anisotropic case we can reconstruct c(x) up to a change of 
coordinates. 

In the classical BC method material parameters are reconstructed using 
hyperbolic techniques, see [H E3, El HB EEJ ES, EZ] . Straightforward numerical 
implementation of the BC method is computationally demanding. Inspired 
by Isaacson's iterative measurement scheme [22] for electrical impedance to- 
mography (EIT), we combine the BC method with an adaptive time reversal 
iteration that lets the medium do most processing of the data. 

Traditional time reversal methods record waves, invert them in time, and 
send them back into the medium. If the recorded signal originated from 
point sources or was reflected from small scatterers in the medium, the time 
reversed waves focus at the source or scatterer points. This is useful for time 
reversal mirrors in communication technologies and medical therapies. For 
early references on time reversal (using ultrasound in air), see the seminal 
works of Fink pjjj CES1 [17]. For a microlocal discussion of time reversal, see 
Bardos (3j[4], and for another mathematical treatment see Klibanov [29J. 

Time reversal in known background medium with random fluctuations 
has also been extensively studied [HI El [7J El UHl E] and applied to medical 
imaging, non-destructive testing and underwater acoustics. These methods 
are outside the scope of this paper. 

Iterative time reversal has been used to find best measurements for inverse 
problems. By "best" we here mean "optimal for detecting the presence of an 
object". This distinguishability problem has been studied for fixed-frequency 
problems in EIT [22] and acoustic scattering [34]. The connection between 
optimal measurements and iterative time-reversal experiments was pointed 
out in [3H [351 EE]. For the wave equation, the best measurement problem 
has been studied in [H], where it has been shown that the optimal incident 
field for probing a half space can be found by an iterative process involving 
time reversal mirrors. 

Let us describe our hybrid method in a simple case. Take a compact set 
M C M. 3 with smooth boundary <9M, and let c(x) be a scalar-valued wave 
speed in M. Consider the wave equation 

u tt - c(x) 2 Au = in M x R+, (1) 
u\ t =o = 0, 
u t \ t =o = 0, 
-c(x)~ 2 d n u = f(x, t) in SMxl + , 

where d n denotes the Euclidean normal derivative. We denote by = 
uf(x,t) the solution of (pp) corresponding to a given boundary source term /. 
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Take T > larger than the diameter of M in travel time metric, so 
that any point inside M can be reached by waves sent from the boundary 
in time less than T. We consider imposing a boundary source / on the set 
dM x (0, 2T) and measuring the boundary value of the resulting wave 
during the time < t < 2T. All such boundary measurements constitute the 
Neumann-to-Dirichlet map 

A 2 T : / | — > U f \dMx{0,2T), 

that models a variety of physical measurements [28] ■ Instead of the full map 
we only need here the values A 2T f for a specific collection of sources / given 
by an iterative process. 

We need three special operators on the function space L 2 (dM x [0, 2T]). 
First, the time reversal operator 

Rf(x,t) = f(x,2T-t), 

second, the time filter 

^ i-m.in(2T-t,t) 

Jf{x,t) = 2 y f(x,s)ds, 

and finally the restriction operator 

Pf(x,t) = XB(x,t)u(x,t), 

where xb is the characteristic function of a set B C dM x [0, 2T}. 
We define the processed time reversal iteration 

F := -P(RA 2T RJ - JA 2T )f, 

a n ■= A 2T (h n ), b n := A 2T (RJh n ), (2) 
a 1 

h n+ i := (1 - -)h n - -{PRb n - PJa n ) + F, 

UJ UJ 

where / G L 2 (dM x [0, 2T]) and a,u > are parameters. Starting with 
ho = the iteration ((2j) converges to a limit h(a) = limn^oo h„. In the case 
B = dM x [T — s, T] we have 

lim u h{a \x, T) = (l - XN(x))u f (x, T) in L 2 (M), (3) 

where N C M is the set of points x G M with travel time distance to dM at 
most s. In other words, N is the maximal region where the waves sent from 
the boundary can propagate in time s. 



3 



Let y G dM and denote by j y>v the geodesic (in travel time metric) 
starting from y in the normal direction and parametrised by arclength. Take 
a source / concentrated near (y, T—s) G dMxM. + and vanishing for t < T—s. 
Then u^(x,T) is supported in the union of the set N and a neighbourhood 
of the point Xo = J y ,v(s), and the limit ([31) vanishes outside a neighbourhood 
of Xo. This way our iteration produces a wave localized near Xq at time T. 

Below we introduce an alternative construction of focusing waves for gen- 
eral sources / by applying the processed time reversal iteration with two 
different sets B and using the difference of the resulting boundary values. 

We show also that the limits (02) can be used to determine the travel time 
distance between an arbitrary point x = j V)V (s) G M and any boundary 
point z G dM. It follows that the processed time reversal iteration gives 
an algorithm for determining the wave speed in the medium. Furthermore, 
given a source /, our iteration allows us to find a wave u h (x, t) having at time 
t = T\ approximately the value Xn(x)u^ (x, T 2 ) for any T 2 without knowing 
the material parameters of the medium. Analogous methods without the 
multiplier Xn(%) have been developed before in [23l l3T] . 

The paper is organised as follows. In Section [2] we formulate our main 
results. In Section [3] we prove the convergence of the processed time reversal 
iteration. Section [4] is devoted to the analysis of the focusing properties of 
the waves, and in Section [5] we show how to reconstruct the metric using the 
boundary distance function. In the last section we discuss our method in the 
case of noisy measurements. 



2. Definitions and main results 

Let us consider the closure M C M m , m > 1, of an open smooth set, or a (non- 
compact or compact) complete Riemannian manifold (M, g) of dimension m 
with a non-empty boundary. For simplicity, we assume that the boundary 
dM is compact. Let u solve the wave equation 

u tt (x,t) + Au(x,t) = in M x R + , (4) 
u\ t =o = 0, u t \ t =o = 0, 

B v ^U\dMxM,+ = f ■ 

Here, / G L 2 (dM x M + ) is a real valued function, A is a formally self-adjoint 
elliptic partial differential operator of the form (in local coordinates in the 
case when M is a manifold) 

m Q ( dv \ 

Av = -J2 /^W^IflWr^ (K x )\g( x )\h jk (%)-^(x)) +q(x)v(x) (5) 
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where [g^ k ] is a smooth real positive definite matrix, cqI < [g^ k {x)} < C\I y 
c<), Ci > 0, \g\ = det([gjk\), where [g^] is the inverse matrix of [g^ k ], and 
ji{x) > C2 > and q(x) are smooth real valued functions. Also, 

B ViV v = —d u v + rjv 

where rj : dM — > R is a smooth bounded function and 

m d 
d u v = ^ K x )9 3k (x)u k —v(x), 

j,k=i ' A ' 

where u(x) = (vi, 1/2, . . . , v m ) is the interior co-normal vector field of dM 
normalised so that Y^k=i9^ kv o Vk = 1- A particular example is the operator 

A = -c 2 (x)A + q(x) (6) 

for which d v v = c(x)~ m+1 d n v, where d n v is the Euclidean normal derivative 
of v. 

We denote the solutions of jl]) by 

w (x, t) = u(x, t). 

For the initial boundary value problem pj we define the response operator 
(or non-stationary Robin-to-Dirichlet map) A by setting 

A/ = U f \ dM xR+- (7) 

We also consider the finite time response operator At corresponding to time 
T > 0, 

A T f = U f \ dM x(0,T)- (8) 

By |33], the map A T : L 2 (dM x (0,T)) -> HV 3 (dM x (0,T)) is bounded, 
where H s (dM x (0,T)) denotes the Sobolev space on dM x (0, T). Below 
we consider A T as a bounded operator that maps L 2 (dM x (0, T)) to itself. 

We also need some other notation. The matrix gjk{x), that is, the inverse 
of the matrix g 3 (x), is a Riemannian metric in M that is called the travel 
time metric. The reason for this is the waves propagate with speed one with 
respect to the metric ds 2 = Yljk 9jk{x)dx : 'dx k . We denote by d(x,y) the 
distance function corresponding to this metric. For the wave equation we 
define the space L 2 (M,dV^) with inner product 

(u,v) L 2 {M4Vfi) = / u(x)v(x)dV fl (x), 

J M 
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where dV^ = (i(x)\g(x)\ 1 / 2 dx 1 dx 2 . . .dx 
For t > and T C dM, let 



M(T,t) = {xeM : d(z,r) < t} 



(9) 



be the domain of influence of T at time t. 

We denote L 2 (B) = {/ G L 2 (<9M x E + ) : supp (/) cB},Bc9Mx E + , 
identifying functions and their zero continuations. When T C dM is open 
set and / G I 2 (r x it is well known (see e.g. [21]) that the wave 

v,f(t) = v,f(- ,t) is supported in the domain M(T,t), 



2.1. Results on convergence of processed time reversal iteration. 

Our objective is to find a boundary source / such that the wave u^(x,T) is 
localized in a neighbourhood of a single point xq. Note that in this paper 
we do not focus the pair (u(x, T), u t (x, T)), but only the value of the wave, 



Let < T < T, T C dM be open set and / G L 2 (T x E + ). Our first aim 
is to construct boundary values h(a) G L 2 (pM x R + ), a > such that 



in L 2 (M). Here Xx(x) is the characteristic function of the set X. Then 



is supported in M(T,T) \ M(dM, T ) (see Figure 1). When T -> T and 
— > {2} then under suitable assumptions discussed later the set M(Tj,T) \ 
M(dM,T ) goes to a single point. Thus for / G L 2 (Tj x the waves 
(1 — XM(8M,T ))uf (T) localize to a small neighbourhood of a single point. For 
later use we will consider (fTUl) below in a more general setting. 

Next we explain a procedure to find boundary values h(a). Denote 



u 



f (t) G L 2 (M(T,t)) = {ve L 2 (M) : supp (v) C M(T,t)}. 



u(x, T). 



]imu h{a >(T) = X M(dM,T )U f (T) 



(10) 



U m u /-ft(«)( T ) = (1 - X M(dM,T ))u f (T) 




(11) 




where J(s,t) = \xL(s,t), 



L = {(s, t) G R + x R + : t + s < 2T, s > t}. 



(12) 
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Figure 1: The set M(T, T) \ M(dM, T ). Our goal is to find waves sent from 
the boundary that focus into this area. 




We call R = R 2 t the time reversal map and J = J 2 t the time filter. We 
denote by P = P B : L 2 (dM x [0, 2T]) -> L 2 (<9M x [0, 2T]) the multiplication 
operator 

Pijf(x,t) = XB(x,t) f(x,t), 

where 5 = U/=i( r i x [ T - Tj,T\), r j C are open sets and < T 3 < T. 
Next, we consider A 2 t, R-it-, and J 2 t as operators from L 2 (dM x [0, 2T]) to 
itself. 

Let a G (0, 1) and w > be a sufficiently large constant. We define 
a n ,b n E L 2 (dMx[0,2T]) and h n = h n (a) E L 2 (dM x [0, 2T]) by the iteration 

a n ■= A 2T (h n ), b n := A 2T (RJh n ), (13) 
a 1 

h n+1 := (1 - -)h n - -(PRb n - PJa n ) + F, 
with a = 0, b = 0, /i = 0, and 

F = -P(RA 2T RJ - JA 2T )f- 

UJ 

We say that this iteration is the processed time reversal iteration on time- 
interval [0, 2T] with projector P and starting point / G L 2 (dM x [0, 2T]). 

Here a n corresponds to the "iterated measurement", b n to the "time filtered 
and time-reversed measurement" and h n+ \ to post-processing of h n , a n and 
b n using time reversal and time filtering. 
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Theorem 1 Let T > 0. Assume we are given dM and the response operator 
A 2 t- Let Tj C dM, j = 1, . . . , J be non-empty open sets, < Tj < T, and 
B = U/=i( r i x [T-Tj'T])- Let f G L 2 (dM x R + ) and let for u large enough 
and a G (0, 1) functions h n = h n (a) be defined by the processed time reversal 
iteration ( 1731) with projector Pb and starting point fo = f\dMx(o,2T)- These 
functions converge in L 2 (dM x 



in L 2 (M), where N = \J J j=1 M(Tj, Tj) C M. 

This theorem is proven later in Section [3l Before that we consider its 
consequences. 

2.2. Results on focusing of the waves Let us consider a geodesic 7^ 
in (M,g) parametrised along arc length where 7^(0) = x, 7r,§(0) = £ with 
||£|| g = 1. Let v = v(z), z G dM be the unit interior normal vector of dM 
in (M, g). There is a critical value r(z) G (0, 00], such that for t < t(z) the 
geodesic 7*,i>([0, £]) is the unique shortest geodesic, and for t > t(z) it is no 
longer a shortest geodesic from its endpoint J z ,v(t) to dM. 

We say that Tj — > {z} if C r\, and f]JL 1 Tj = {z}. 

Theorem [T] yields the following results telling that we can produce focusing 
waves, that is, wave that focus to a single point. 

Corollary 2 Let z G dM and < T < f < T . Let x = yg,„(T) and 
Tj C dM, j G Z + be open neighbourhoods of z G dM such that Tj — > {z} 
when j — > 00. 

Let f G C^°(dM x IR+). Let h n (a;T ,j) be the functions obtained from 
the processed time reversal iteration (Q3|) with projector Pb, B = (Tj x [T — 
T,T]) U (dM x [T — To,T]) and starting point fo = /|aMx(o,2T)- Similarly, 
let h' n (a;T 0l j) be the functions obtained from the processed time reversal 
iteration (Q5]J with projector Pb< , B' = dM x [T — T ,T] and starting point 



h(a) = lim h n (a) 



and the limits satisfy 



\imu h{a) (x,T) = XN (x)u f (x } T) 



fo- Let 



h n (a;T ,j) = f - h n (a;T ,j) + h' n (a;T Q , j). 



IfT < t(z) then 



(14) 



lim lim lim lim 

T ^T a ~*0 n—*oo 



1 



(^)(T) = C (x)uf(x,T)5^(x) 



u 



(T - T )( m+1 )/ 2 
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Figure 2: Subsets of M used in the proof of Corollary [3 

x 



















T 



dM 2 Tj 



in T>'(M) where C (x) > does not depend on f. IfT> r(z) the limit ( [7^ ) 
is zero. 

In Figure^ the set M(Vj,f) \ M(dM, T ) tends to the point x when 
j — » oo and T — > T. This turns out to be essential in the proof of Corollary 

El 

Above, 5% is Dirac's delta distribution on (M, g) such that 
Jm 

2.3. Inner products. Later, we show that that using f,h G L 2 (dM x 
[0, 2T]) and the boundary measurements A 2 t we may compute via an explicit 
integral the inner product 

/ u f {x,T)u h {x,T) dV^x) = I {Kf){x,t)h{x,t)dS g {x)dt, (15) 

JM JdMx[0,2T] 

where dS g is the Riemannian surface volume of dM and 

K = Kit '■= R.2T^-2tR2tJ2T ~ J2T^-2T (16) 

so that terms (PT5l) can be found using measurement operator A 2 t and simple 
basic operations like time reversal i?2T and the time filter operator J^t- 

In the above formula the Riemannian surface volume of dM can be com- 
puted using the intrinsic metric of dM , which is determined by the map A 2 t- 
Indeed, its follows from Tataru's unique continuation principle (see [381 [40] . 
see also Theorem [4] below) that the Schwartz kernel of A 2 t is supported in 
the set E = {(x,t,x',t f ) E (dM x [0,2T]) 2 : t - t' > d{x,x')} and that 
the boundary dE is in the support. The set dE determines the distances 
of points z,z' e dM with respect to the metric of (M,g), and thus also the 
distance with respect to the intrinsic metric of the boundary (dM,ggM)- 
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2.4. Reconstruction of material parameter functions. It is well known 
that the response map A can not generally determine coefficients fj, and g^ 
because of two transformations discussed below. 

First, we one can introduce a coordinate transformation, that is a diffeo- 
morphism F : M — > M such that the boundary value F\qm is the identity 
operator. Then the push forward metric g = F*g, that is, 

m(y) = f^f^ fe( x )> y = F ( x ) 

p,q=l ^ ^ 

and the functions Jl = /ioF" 1 , q = qoF' 1 , and rj = r] determine the operator 
A of the form ((5j) such that response operators for A and A are the same (for 
this, see [251 M, E7J). 

Second, one can do the gauge transformation u(x,t) — > K,(x)u(x,t) where 
K G C°°(M) is strictly positive function. In this transformation the operator 
A is transforms to the operator A K defined be 

A K w := kA(— w) 

K 

and the operator in boundary condition is transformed to B u ^, with 

rj = rj — n~ 1 d u K. 

Then the response operator A K of A K coincides with the response operator A 
of A 

By above transformations, making a gauge transformation with k = fi , 
we come to operator A K having the same response operator as A and that 
can be represented in the form ((HJ) with \i{x) — 1. It turns out that this is 
the only source of non-uniqueness. 

Corollary 3 Assume that — \ and that we are given the boundary dM and 
the response operator A. Then using the the processed time reversal iteration 
we can find constructively the manifold (M, g) upto an isometry and on it 
the operator A uniquely. 

Moreover, if M C M. m , m > 2, and A is of the form (EJJ, given the set M 
and the response operator A we can determine c(x) and q{x) uniquely in M . 

We note that the unique determination of (M, g) and A is well known, see 
e.g. [HI EHJ [26] and references therein. The novelty of Corollary [3] is that in 
the construction based on processed time reversal iteration the use of iterated 
measurement avoid many computationally demanding steps used required in 
[26] . An example of such step is a Gram-Schmidt orthogonalisation of basis of 
spaces L 2 (r x [t 1 ,t 2 ]), T C dM with respect to an inner product determined 
by A 2 r- This is a typical step used in traditional Boundary Control method 
that is very sensitive for measurement errors. 
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3. Proof of the convergence of the iteration 

3.1. Controllability results. The seminal result implying controllability 
is Tataru's unique continuation result (see [381 SO]) 

Theorem 4 (Tataru) Let u be a solution of wave equation 

u t t{x, t) + Au(x, t) = 0. 

Assume that 

w|rx(o,2r) = 0, d u u\ rx{0t 2r) = (17) 

where T C dM, F ^ is open, and r > 0. Then, 

u(x, t) = 0, d t u(x, r) = for x G M(T, r). 

This result yields the following Tataru's controllability result, see e.g. [26] 
and references therein. 

Theorem 5 Let T C dM be open and r > 0. Then the linear subspace 

{t*'(r) G L 2 (M(r,r)) : / G C °°(r x [0,r])} 
zs dense in L 2 (M(r,r)). 

3.2. Inner products. Let us recall the Blagovestchenskii identity: 
Lemma 6 Let f,he L 2 (dM x [0, 2T]). Then 

[ u f (x,T)u h (x,T)dV fM (x) = (18) 

J M 

[ [ J(t, s) [f(t)(A 2T h)(s) - {A 2T f)(t)h(s)] dS g (x)dtds, 

J[0,2T] 2 J8M 

where J(t,s) = \xi(s,t), see fj^). 

The proof in slightly different context is given e.g. in |26j. 

Proof. Let w(t, s) = f M u*{t)u (s) dV^. Integrating by parts, we see that 

{dt-dl)w{t,s) 

= - [ [Au f {t)u h {s) - u f {t)Au h {s)] dV^x) 

= - [ [d„u f {t)u h {s) - u f {t)d u u h {s)] dS g (19) 

JdM 

= [ [(-d v uf(t) + V u f (t))u h (s) - u f (t)(-d u u h (s) + vu h (s))] dS g 

JdM 

= [ [f(t)A 2T h(s) - A 2T f(t)h(s)] dS g . 

JdM 
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H=o = W \ S=Q = 0, d t w\ t=0 = d s w\ s=0 = 0. 



Moreover, 



Thus we can consider ( 11911 as one dimensional wave equation with known 
right hand side and vanishing initial and boundary data. Solving this initial 
boundary value problem we see the claim. 

□ 



Consider now the operator K appearing in formula ( 1151 ) in more detail. 
The Schwartz kernel of A 2 t is the Dirichlet boundary value of the Green's 
function G(x, x', t — t') satisfying 

(d 2 t + A)G x , tV (x,t) = 6 x >(x)6(t-t') inMxM + , (20) 
G x ' t t'\t=o — 0, d t G x i tt '\t=Q = 0, 
B u ^G x ' :t '\dMxR + — 0, 

where G x i >t i(x, t) = G(x, x', t — £'). Then 

G(x, x', t-t')= G(x', x, t - t') 



and we see that 



Any — -R2TA2T-R2T 



where R2rf(x, t) = f(x, 2T — t) is the time reversal map. 

Thus, as A-2 T = R 2 t^-2tR2T^ we can write result of Lemma [6] as 

u f {x,T)u h (x,T)dV t t(x) = [ (Kf)(x,t)h(x,t)dS g (x)dt (21) 

M JdMx[0,2T] 

where K is defined in (TTBl). 



3.3. Proof of Theorem [U. To produce focusing waves, let us consider the 
problem 

h ^ B) h f (T)-u\T)\\l 2{MjdVii) (22) 

where B = \J J j=1 {Tj x [T - T h T}), Tj C dM are open and < Tj < T. 

As the solution of this minimisation problem does not always exist, we 
can regularise the problem and study 

min F(h,a) (23) 

h£L2(dMx[0,2T]) 
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where a G (0, 1) and 

F(h, a) = (K(Ph - f), Ph - f) L * {dM x[o,2T],ds g ) + oi\\h\\ 2 L 2 {dMx[0t2ndSg) . 

We recall that P = Pb is multiplication with the characteristic function of 
P, that is, {PbH){x, t) = Xb(x, t) ■ h(x, t). 
By (EH), 

F(h,a) = \\uf(T) -u ph {T)f L2{M4Vii) + a||/i||! 2(9Mx[0i2T]) . 
The minimisation problem (l23l ) is equivalent to (l22l) when a = 0. 

Lemma 7 For given a G (0, 1) the problem (d3Jj has a unique minimiser. 
Moreover, the minimiser is the unique solution of the equation 

(PKP + a)h = PKf. (24) 

Proof. The minimisation problem is strictly convex, and the map h i— > u h (T) 
is continuous L 2 {dM x [0,2T]) -> Pf Sl (M), Sl < 5/6 by [33j, see also [H2]. 
Now if hj -> ft weakly in L 2 (<9M x [0,2T]) as j oo then Pftj -> Pft 
weakly in L 2 (dM x [0, 2T]). As the embedding H Sl (M) -> L 2 (M) is compact 
and the map ft i-> it ft (T) is continuous L 2 (dM x [0, 2T]) -> Pf Sl (M), then 
u^i(T) -> m p,1 (T) in L 2 (M). Therefore, 

lim {K(Phj - /), Pftj - / Wear* [0,221) 

Let ftj be a sequence such that lim. ? _ ) . 00 F(hj, a) = inf/ l P(ft,a). As 
F(h,a) > 

& II ft II \i i the sequence (ftj) is bounded. Since ||ft||^2 < liminfj^oo ||ftj||^2, 
this implies that F(h,a) < liminf F(hj, a). 

Thus by choosing a subsequence of ftj, we can assume that hj converge 
to ft in the weak topology of L 2 [pM x [0,2P]), and the limit ft is a global 
minimiser of P(- , a). 

By computing the Frechet derivative of the functional P(ft, a) in any 
direction ft G L 2 [dM x [0, 2T]) at a minimiser ft, we see that 

= D h ((K(Ph - f), Pft - /} + a||ft|| 2 2 )ft 
= (ft, P*K{Ph - /)) + (ft, P*K*{Ph - /)) + 2a(ft, ft) 

for all ft. Since if* = if and P* = P, this implies that ft satisfies the equation 
(1241 ). As If is non-negative, PKP + al > al so that equation ( 1241) has a 
unique solution providing the minimiser ft = ft(a). □ 
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Note that equation (1241 also implies that 

h(a) G Ran(P). (25) 



Next we want to solve equation (1241 ) using iteration. To this end, let 
uj G M+ be a constant such that 

uj > 2(1+ \\PKP\\). 

Then the equation (1241) can be written as 



1 a + PKP 
(/ - S)h = - PKf, where S = I- ^—^ . (26) 

Since PKP is non-negative and 

al<al + PKP < -I, 
- 2 ' 

we see that \\S\\ < 1. Thus we can solve h using iteration: Let 

F := -PKf = -P{JA 2T ~ RA 2 TRJ)f, 

to CO 

ho = 0, and consider the iteration 

h n+1 = Sh n + F, n = 0, 1, — 

As h n = Ph n by (1251) . we can write the iteration in the form ( Tl"3l) . and 
lim^oo h n = h(a) in L 2 {dM x [0, 2T]). 

Applying this algorithm can find the minimisers h = h(a) G Ran (P) of 
problem (1241) . They converge by the following lemma: 



Lemma 8 We have 

\\mu h{a \x,T) = XN(x)u f (x,T) 

in L 2 (M), with N as in TheoremUl 
Proof. Clearly, 

F(h,a) = F + \\x N (u Ph {T) - u f (T))\\ 2 L3{M) + a||/i||i2 (9Mx[0)2T]) , 

^o=||(l-X^ / (T)||i 2(M) . (27) 

Here the first term does not depend on h. By Theorem El for any e > there 
is a h e e L 2 (B) such that 

\\u P HT)-XNU f (T)\\l 2(M) <^. (28) 
14 



Thus we have 



e 



F(h e ,a) = F + - + a||/i(e)||i2 (aMx[0)2T]) . 



This shows that if a < a(e), where 



a e 



L 2 (dMx[0,2T}) 

then 

F(h e , a)<F + e. 

Thus the minimiser h(a) of F(h, a) with a < a(e) satisfies 

F(h(a),a) <F + e. 

As by ((251), h(a) = Ph(a) so that h(a) G L 2 (r x [0,T]), we get from (l28il 
that, for a < ct(e), 

||^(^(T)-^(T))||i 2(M) <e/2. (29) 

As supp (w h ( Q )(T)) C N, the claim follows. □ 

Theorem Q] follows from Lemmata [7| and El □ 

4. Proofs for the focusing of the waves 

Next we prove Corollary [2l 

Proof of Corollary [U. Let {z(x), s(x)) be the boundary normal coordinates 
of x G M, that is, s(x) = d(x, dM) and z(x) is the closest point of DM to x 
when such a point is unique. When a closest boundary point is not unique, 
the boundary normal coordinates are not defined. 

We consider first the claim of the corollary in the case when T < t(z). 
Then the boundary normal coordinates near x = 72^ (T) are well defined and 
the metric tensor in these coordinates has the form 

9=1 I r„ J, oM«»-i J > ( 30 ) 



[<M^)]™Si 

where [(^(z, s)] G R(" 1 ~ 1 ) x ( m_1 ) is a smooth positive definite matrix valued 
function near (z,T) = (z(x),s(x)). Then there is a Ci = Ci(z,T) > such 
that 

Cf 1 /^ [fctfCz,*)] <CiJ 
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near (z, T). Using this we see that there is a C 2 > such that 

c _! < vol g (M({*},f)\M(0M,r o )) < c 
2 ~ (f -T )(f 2 -T 2 )( m - 1 )/ 2 ~ 2 ' 

As the solution u f (x,t) is smooth for / e C£°(dM x R + ), 

^ XM( {?} ,f)UM(OM,T ) - XM( 9 M,T ) ^ = ^ 

To^f vol,(M({z},T) \ M(«9M,T )) 

in V'(M), the claim follows in the case T < r(z). 

When T > t(z), the claim follows from the fact that for sufficiently small 
f - T > the set M({z},f) \ M(dM, T ) is empty. □ 

Remark The above proof of Corollary [2] also shows that in Corollary [2] 

(f - T' n V m + 1 )/ 2 

C (x ) = lim — ^ -. (31) 

To ^f vol, (M({z},T)\M(dM, T )) 

5. Boundary distance functions and reconstruc- 
tion of the metric 

Here we present a method for finding the boundary distance functions using 
processed time reversal iteration. 

Let z,y G dM, < T\ < t(z) and T > dg M (y,z) + T\. Denote x = 
lz,v(Ti). Next we give an algorithm that can be used to determine d(x, z). 

To this end, let Tj C dM and Ej C dM be neighbourhoods of z and y, 
respectively, such that IT, — > {z} and Ej — >■ {y} when j — >■ oo. 

Let e > 0, r e [0,T], and 

N] = M(Tj, TO B] = r j x [T - Tx, T] 

iVf = M(E i; r) B? = Ej x [T - r, T] 

iV 3 = Af (9M, Ti - e) B 3 = <9M x [T — (7\ — e), T]. 



Lemma 9 The distance d(x, z) is the infimum of all r e [0, T] i/iai satisfy 
the condition 

the set J(j, e) = (Nj n iV 2 ) \ A" 3 contains (32) 
a non-empty open set for all j G Z+, e > 0. 
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See Figure [3] for a sketch of the cases of empty and non-empty I(j,e). 

Proof. First consider what happens if d(x,y) < r. Since Ti < r(z), we see 
that then B(x,r) fl (Nj \ A^ 3 ) contains a non-empty open set for all r > 0, 
where B(x,r) C M is a ball of (M,g) with center x and radius r. When 
r < r — d(x,y), we see that B(x,r) C iVj. Thus I(j,e) contains an open set 
and (l32l) is satisfied. 

On other hand, if d(x,y) > r, let r = d(x,y) — r. When j — > oo and 
e — > 0, we see using metric in the boundary normal coordinates (l30l) that 
A^ 1 \ A^ 3 — > {x} in the Hausdorff metric. Thus when j is large enough and e 
is small enough, Nj \ c B(x, r/2). Then B(x, r/2) n Nf = 0, and {32j) is 
not satisfied. 

Summarizing, condition (|32l) is satisfied if <i(a;, y) < r and not satisfied if 
d(x,y) > r. As d(x,y) < T, this yields the claim. □ 

Next we show that by using the processed time reversal iteration, we can 
test if condition (l32l ) is valid. To this end, denote 

N, (j, 6) = N}UN* B 1 (j,e) = B) U B\ 

iV 2 (j, e) = Nj U Nf B 2 (j, e) = B) U B\ 

N 3 (j, e) = NjuNfu N* B 3 (j, e) = B) U B? U B\ 

N^e) = N? B 4 (j,e) = B 3 e . 

Let / G C^°(dM x R+). Using the processed time reversed iteration on 
time interval [0, 2T] with projectors Pb corresponding to B = Bk(j,e), k = 
1,2,3,4 and starting point /, we obtain functions h n (a;e,j,k) G L 2 (dM x 
[0,2T]). Using them, define 

p n (a,j,e) = h n (a;e,j, 1) + h n (a;e,j,2) - h n (a;e,j, 3) - h n (a; e,j, 4). (33) 
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Lemma 10 The condition (CHjj is satisfied if and only if there exists an 
f G C^(dM x R + ) such that for any j G Z + and e > the functions 
p n (a,j,e) defined in formula (SB) satisfy 

lim lim (K 2T p n (a,j,e),p n (a,j,e)) ^ 0. (34) 

a— >0 n— +oo 

Proof. The functions h n (a; e, j, k) defined by processed time reversal itera- 
tion satisfy 

XN k(j e) uf ( T ) = lim n lim u h ^> k \T), k = 1,2,3,4. 
A simple computation gives us 

XlU,e)( X ) = XNx(j,e)( X ) + XN 2 (j,e)( X ) ~ XN 3 (j,e)( X ) ~ XN 4 { h e)( X ) ( 35 ) 

for all x G M. Therefore, using (l33l we see that in L 2 (M) 
Xi<je)U f (T) = lim lim u Pn{a ' j ' e) (T). 

By Theorem [5] we see that the functions u^(T) with / G C£°(dM x R + ) are 
smooth functions that form a dense set in L 2 (M(<9M, T)). Thus condition 
(l32l) is satisfied if and only if there exists an / G C^(dM x R + ) such that 
for any j G Z + and e > 

(Xl{j,e)U f (T),U f (T)) L 2 {M) 

= lim lim {K 2T p n {a,j, e),p n (a,j, e)) L 2r dMx[02T n ^ 0. 

This proves the claim. □ 

We note that if condition (1341) is satisfied for some / G C£°(dM xR + ), 
it is satisfied for all such / in an open and dense set. 

Lemmata l9l and [POl give an algorithm for the determination of d(y, x) by 
using processed time reversal iteration 

d(y, x) = inf{r G [0, T] : there is an / G C™{dM x R + ) (36) 
such that (pm holds for all j G Z + and e > 0}. 

Summarizing, we have proven: 

Theorem 11 Assume we are given dM and the response operator A. Let 
z,y G dM, T\ < t(z). Then using the algorithm [3b)) we can compute d(x, y) 
for x = 7 Z)1 /(Ti). 
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Let us consider consequences of the result above. To this end, we define 
the set of the boundary distance functions. For each x G M, the correspond- 
ing boundary distance function, r x G C(dM) is given by 

r x :dM-> R + , r x (z) = d(x, z), z G dM. 

In fact, r x G Lip (dM) with the Lipschitz constant equal to one. The bound- 
ary distance functions define the boundary distance map 1Z : M — > C(dM), 
TZ(x) = r x , which is continuous and injective (see [26J ) . Denote by 

U(M) = {r x G C(dM) : x G M}, 

the image of K. It is known (see [261 [27]) that, given the set TZ(M) C C(dM) 
we can endow it, in a constructive way, with a differentiable structure and a 
metric tensor g, so that (TZ(M),g) becomes a manifold that is isometric to 

(M,g), 

(1Z(M),g) = (M,g). 

See [27J. A stable construction procedure of (M,g) as a metric space from 
the set TZ(M),g) and Holder type stability estimates are give in [24|. 

Example By the triangle inequality, 

\\r x -r y \\ 00 <d(x,y). (37) 

We consider in this example the case when (M, g) is a compact manifold 
such that all points x,y G M can be joined with a unique shortest geodesic. 
This implies that, for any x, y G M, the shortest geodesic 7([0, s}) from x 
to y, parametrised along arc length, can be continued to a maximal geodesic 
7QO, L}) that hits the boundary at a point z = j(L) G dM. Then 

\ r x(z) - r y (z)\ = d(x,y) 

implying equality in (1371) . Thus in the case when all geodesies between ar- 
bitrary points x,y G M are unique, the manifold (M, g) is isometric to the 
manifold 1Z(M) with the distance function inherited from C(dM). In the 
general case the construction of the metric is more elaborate. 

By Theorem [TT] we can compute for all x = j ZtV (T) with T < r(z) the 
corresponding boundary distance function r x . Since all points x G M can be 
represented in this form (see e.g. [15]) we can find the set 1Z(M) that can be 
endowed with a manifold structure isometric to the original manifold (M,g). 
We have thus proven the following result: 
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Corollary 12 Assume we are given dM and the response operator A. Then 
using the processed time reversal iteration we can find the manifold (M, g) 
upto an isometry. 

Corollary [12] can also be formulated by saying that we can find the metric 
tensor g in local coordinates. For example, for any point Xq G M there 
are Zj G dM, j = 1,2,..., m, such that x t— > (X 1 (x), . . . , X m (x)) with 
Xi(x) := d(x,Zj) define local coordinates near xq. In these coordinates the 
distance functions x \— > d(x,z), z G dM determine the metric tensor. For 
details of this construction, see [26] . 

Next we prove Corollary [3l 

Proof of Corollary [3]. As the boundary data dM and A determine (M, g) 
upto isometry, we can apply formula (|3TT) to find the function C (x), x G M 
in Corollary [2] in any local coordinates of (M,g). Thus in local coordinates 
we can find the values of waves u'{x,t) for all x G M, t > 0. By Tataru's 
controllability theorem, the waves u^(x,T) with fixed T > form a dense 
set in L 2 (M(dM,T)). In this dense set we can find in local coordinates the 
distributions 

Au f (x,t) = -d$u f (x,t) = -u fu (x,t), 

implying that we can find values Aw for all w G L 2 (M(dM,T)). As T is 
arbitrary, we find Aw for all w G L 2 (M). 

The case when M C M. m , m > 2 and A is of the form (jOj), we use the fact 
that a manifold with conformally Euclidian metric is embedded into M. m in 
a unique way when dM is fixed. Thus we find c(x) when a metric g is given, 
and as above we also find an operator A K that is equivalent to A through a 
gauge transformation. As there is a unique operator of the form ((6]) in the 
orbit of the the gauge transformations A —> A K the claim follows (see [26]). 
□ 



6. Iteration when measurements have errors 

Let (O, E, P) be a complete probability space. 

Assume the measurements have a random noise, that is, measurements 
give us with an input / the function A 2 rf + e, where e is random Gaussian 
noise that has values in L 2 (dM x (0, T)). Assume that Ee = and denote 
the covariance operator of e by C f . Note that C € is a compact operator on 
L 2 (dM x (0, T)) (see [12]) and thus the standard white noise on dM x (0, T) 
does not satisfy our assumptions. 
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Assume that the noise is independent of previous measurements each time 
when we do a new measurement. When the noise is added to the processed 
time reversal iteration, we come to the iteration of the form 

h n+1 = Sh n +p + N n , 

where N n = PJe^ — PRe\ where e l n and e\ are random variables having the 
same distribution as e. Thus the N n are independent identically distributed 
Gaussian random variables satisfying KN n = and having covariance oper- 
ator C N = PJC e J*P* + PRC e R*P*. Let us consider the averaged results of 
iterations 

i K „ 

K e = x E ( 38 ) 



Then 

K , K „ K 



- sy 2 N n 



K 

n 




n+2 - {n + 2){I - SY 1 S n+1 )N n 



H K + H K + H Ki 



where h n are the results of the processed time reversal iteration without noise. 

Above, the deterministic term H\ converges to the same limit lim^oo h n = 
h{a) as the processed time reversal iteration without noise. Now consider PP\ 
and H\ as random variables in L 2 (dM x (0, T)). They can also be viewed 
as random fields on dM x (0, T), see [37]. By the law of large numbers in 
infinite dimensional spaces [20], we see that 



lim \\H 2 K \\ L 2t aMx r^ T) \ = 0. (39) 

As 1 1 jS" 1 1 x,(x, 2 (aJW)) ^ \, the last term H\ satisfies also an estimate analogous 
to (l39l) . Thus the averaged processed time reversal iteration with noise con- 
verges to the same limit as processed time reversal iteration without noise, 
that is, 

lim h°£ e = h{a) in L\Q; L 2 (dM x (0, T))). 

K^oo 
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